f/s  i*/i 


.  AD-A110  5*3  TEXAS  TECH  UNIV  LUBBOCK  DEPT  OF  MATHEMATICS 

APPLICATIONS  OF  INVARIANT  IMSEDOINS.(U) 

DEC  SI  p  NELSON  AFOSA-B0-0aSS 

UNCLASSIFIED  AF0SR-TA-B2-002B  NL 


Kea>n  ADAH  054  3 


AP0SR-m-  a2  - 


002  8 


LEVEL* 


d^- 


ir 


r  r.  a  iv, 


*!. ""  •  *  r*n  t-,  ; 

v';  f  Lo  J  •«* 


FINAL  SCIENTIFIC  REPORT 

Air  Force  Office  of  Scientific  Research  Grant  AFOSR-80-0235 
Period:  1  July  1980-30  June  1981 
Title  of  Research:  Applications  of  Invariant  Imbedding 
Principal  Investigator:  Paul  Nelson 


49 


Approved  for  public  reloa30  { 
distribution  unlimited. 


UNCl.ASSit  jiID 

m  - - - - - - - - - - — 

SECURITY  CLASSIFICATION  OF  THIS  PAGE  Clt'.-.r/,  t>„ i„  Fi.lftudi 


REPORT  DOCUMENTATION  PAGE 


I.  REPORT  NUMBER  /\  ft  O 

AFOSR-TR-  8  2  -  0  0  2  8 


3  A  pc  Ki.AU  IXSTRl'CTMXs 

MOl~ _ bi'.i  oKK  rov;u.i:;,N'i  ;•  okm 

2.  GOVT  ACC  ESSON  NO.'  2  RECIPIEnT-S  CATALOG  NUMBER 


il  6  svz 


14.  title  (and  Subtitle) 


APPLICATIONS  OF  INVARIANT  IMBEDDING 


|7.  AOTHORfs) 


5.  TYPE  REPORT  &  PERiOO  COVERED 

FINAL,  1  JUL  80  -  30  JUN  81 


6.  performing  org.  report  number 


8.  CONTRACT  or  GRANT  NCMuERfs; 


Paul  Nelson1 


AF0SR-80-0235 


9.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS  10.  PROGRAM  ELEMENT.  PROJECT,  T  ASA 

_  .  .  „  ,,  ,,  ,  .  AREA  6  WORK  UNIT  NUMBERS 

Department  of  Mathematics 
Texas  Tech  University 

Box  4139.  Lubbock  TX  79409 _ 61102F;  2304/A3 _ 

II.  CONTROLLING  OFFICE  NAME  AND  ADDRESS  12.  REPORT  DATE 

Mathematical  &  Information  Sciences  Directorate  DEC  81 

Air  Force  Office  of  Scientific  Research  '3-  number  of  pages 

Bolling  AFB  DC  20332 _ 16 _ 

14.  MONITORING  AGENCY  NAME  &  ADDRESS’,'//  different  from  Controlling  Olhce)  15-  SECURITY  CLASS,  (ot  this  report) 

UNCLASSIFIED 

1 5a  DECL  ASSI  F|C  ATI  ON  DD~WN  GRADING” 
SCHEDULE 

16.  DISTRIBUTION  STATEMENT  (o I  mis  Report)  ~ 

Approved  for  public  release:  distribution  unlimited. 


'7.  DISTRIBUTION  STATEMENT  (of  the  abstract  entered  in  Dio ck  20,  i!  different  from  Report) 


ia.  supplementary  notes 


f  19.  KEY  WORDS  (Continue  on  reverse  side  if  necessary  and  identify  by  block  number) 


%(,  ABSTRACT  (Continuo  on  rovers o  side  If  necessary  onJ  identify  by  block  nt.mbvr) 

'Invariant  imbedding  is  a  method  for  the  computational  solution  of  tv;o-point 
boundary-value  problems.  One  common  source  of  such  problems  is  in  the  applica¬ 
tion  of  the  method  of  lines,  or  an  expansion  procedure,  to  a  system  of  partial 
differential  equations.  Such  partial  differential  equations  might  describe, 
for  example,  the  dynamic  mechanics  of  structures  such  as  an  aircraft  fuselage 

or  a  missile  silo.  _ 

(CONTINUED) 


1473  EOlTlON  OF  I  NOV  65  is  OBSOLETE 


JilLJLCc 


_ UNCLASSIFIED _ 

SECURITY  CLASSIFICATION  OF  This  PAGE  f.OxYl  Do;*  Enttfd) 


UNCLASSIFIED 


secup.t'.  ci 


j- ;  |  ,ri  I !  ,1  /  /.  .*»'  (■ 


ITEM  #20,  CONT.: 

ff 

,  .^Phe  research  conducted  under  this  grant  has  been  directed  toward  two  objectives, 
/  The  first  objective  was  to  improve  the  effectiveness  and  efficiency  of  invariant 
imbedding  by  providing  means  for  automatically  controlling  the  associated  compu¬ 
tational  effort.  The  second  objective  was  to  extend  th/fe  range  of  applicability 
of  invariant  imbedding  to  include  singular  two-point  boundary-value  problems. 
Such  singular  problems  arise,  for  example,  by  applying  the  method  of  lines  to 
partial  differential  equations  in  spherical  or  cylindrical  coordinate  systems. 


In  regard  to  the  first  of  these  objectives,  a  number  of  different  methods, 
which  are  generically  termed  'relative-error  monitors',  have  been  developed  and 
are  in  the  process  of  being  subjected  to  computational  experimentation.  As 
regards  the  second  objective,  it  has  been  shown  ho-.'  to  apply  the  method  of 
invariant  imbedding  to  'homogeneous'  linear  two-point  boundary-value  problems. 
The  further  extension  to  inhomogeneous  problems  is  presently  being  pursued.  3 

.  i 

1 


1 


UNCLASSIFIED 


CONTENTS 


Abstract 

Research  Objectives  and  Status 
Automatic  Error  Control 
Singular  Problems 
References 

Publications 

Professional  Personnel 

Interactions 


AIR  FORCE  C3T-rl- 
NOTICE  OF  V7'  ';  r  •••  •'  - 

This  tech:- 

approved  l  _ 

Distrifcv.ticv' 

VuiT'iiliTA  J. 

Chief,  Technical  iPlorr-Ki 


»T pTC  FEsE'.iRCH  fi?SC) 

■-  ^  ■* 

^  \  <?’v°d  ?nd  is 

2  j/. 

i 

i  • 

Mon  Division 


ABSTRACT 


Invariant  imbedding  is  a  method  for  the  computational  solution 
of  two-point  boundary-value  problems.  One  common  source  of  such 
problems  is  in  the  application  of  the  method  of  lines,  or  an  ex¬ 
pansion  procedure,  to  a  system  of  partial  differential  equations. 
Such  partial  differential  equations  might  describe,  for  example, 
the  dynamic  mechanics  of  structures  such  as  an  aircraft  fuselage 
or  a  missile  silo. 

The  research  conducted  under  this  grant  has  been  directed  to¬ 
ward  two  objectives.  The  first  objective  was  to  improve  the  ef¬ 
fectiveness  and  efficiency  of  invariant  imbedding  by  providing 
means  for  automatically  controlling  the  associated  computational 
error  within  a  minimal,  or  at  least  reasonable,  degree  of  computa¬ 
tional  effort.  The  second  objective  was  to  extend  the  range  of 
applicability  of  invariant  imbedding  to  include  singular  two-point 
boundary-value  problems.  Such  singular  problems  arise,  for  example, 
by  applying  the  method  of  lines  to  partial  differential  equations 
in  spherical  or  cylindrical  coordinate  systems. 

In  regard  to  the  first  of  these  objectives,  a  number  of  dif¬ 
ferent  methods,  which  are  generically  termed  relative- error  moni¬ 
tors,  have  been  developed  and  are  in  the  process  of  being  subject¬ 
ed  to  computational  experimentation.  As  regards  the  second  objec¬ 
tive,  it  has  been  shown  how  to  apply  the  method  of  invariant  im¬ 
bedding  to  homogeneous  linear  two-point  boundary- value  problems. 

The  further  extension  to  inhomogeneous  problems  is  presently  being 
pursued. 


ii 


RESEARCH  OBJECTIVES  AND  STATUS 


The  work  described  here  was  concerned  with  automatic  error 
control  for  invariant  imbedding  and  with  the  extension  of  invariant 
imbedding  to  permit  solution  of  singular  linear  two-point  boundary- 
value  problems.  Progress  on  these  two  topics  is  discussed  separa- 
ely  in  the  following  two  subsections. 

Automatic  Error  Control 

The  solutions  of  the  initial-value  problem  which  typically  is 
integrated  numerically  in  applications  of  invariant  imbedding  are 
usually  not  themselves  the  object  of  ultimate  interest.  Rather, 
these  are  intermediate  quantities  which  are  used  via  linear  alge¬ 
braic  computations  to  approximate  the  solution  of  an  underlying 
linear  two-point  boundary-value  problem.  Because  of  this  situation, 
it  is  not  readily  apparent  how  to  use  modern  adaptive  initial-value 
integration  techniques  efficiently  and  effectively  within  invariant 
imbedding. 

Nelson  and  Wiggins  [1]  studied  one  method  of  error  control  for 
initial-value  integrations  which  was  specifically  designed  for  the 
version  of  invariant  imbedding  due  to  M.R.  Scott  [2-5] .  In  the 
work  reported  here  we  view  this  earlier  approach  as  one  of  a  class 
of  methods,  each  of  which  is  defined  by  quantities  which  we  term 
relative-error  monitors  and  associated  sensitivity  functions .  Each 
such  set  of  quantities  provides  a  mechanism  for  approximately  con¬ 
trolling  that  part  of  the  error  in  the  approximate  solution  of  the 
underlying  two-point  boundary-value  problem  which  stems  from  the 
initial-value  integration.  We  now  describe  this  approach  to  error 


2. 


control  in  some  detail. 

Consider  a  linear  two-point  boundary-value  problem  (TPBVP)  of 
the  form 

u'  (z)  =  A  ( z )  u  ( z )  +  B(z)v(z)  +  s+(z),  (la) 

v*  (z)  =  C(z)u(z)  +  D(z)v(z)  +  s_(z),  (lb) 

u (0)  =  a,  v (x)  =  b.  (lc) 

Here  u  and  v  are  respectively  m  and  n-  dimensional  vectors  of  de¬ 
pendent  variables,  A,B,C  and  D  are  given  piecewise  continuous  mat- 
tices  of  the  appropriate  dimensions,  a  and  b  are  given  vectors,  and 

x>0.  If  0  <  x  <  x  ,  where  x  is  the  smallest  value  of  x  such 
—  c  c 

that  (1)  fails  to  have  a  unique  solution  (i.e.  xc  is  the  smallest 
eigenlength  of  (1) ) ,  then  the  solution  of  (1)  is  given  by 

v(z)  =  T_1(z){T(x)b  +  q^x)  -  q^z)},  (2a) 

u(z)  =  R(z)v(z)  +  qr(z)/  (2b) 


where  R,T,qr  and  q^  are  determined  by  the  initial-value  problem 
(IVP) 


R' 

-  B  +  AR-RD- 

■RCR, 

R(0)  =  0, 

(3a) 

T' 

«  -T (CR+D) , 

T  (0) 

=  I  (=  identity  matrix) , 

(3b) 

1 

qr 

=  (A-RC)qr 

-  Rs_ 

+  s+/  qr(0)  *  a, 

(3c) 

*1 

=  -T(Cqr  + 

s_)  , 

qa(0)  =  0. 

(3d) 

In  applying  (2)  to  the  computational  solution  of  the  TPBVP  (1) , 
one  normally  integrates  the  IVP  (3)  from  z=0  to  z=x,  with  the  values 
of  R,T,qr  and  qr  at  reporting  points  (=  values  of  z  at  which  it  is 

ultimately  desired  to  compute  u  and  v)  being  stored  during  this 
integration.  Once  T(x)  and  q^ (x)  have  been  computed,  the 


approximations  to  u  and  v  at  the  reporting  points  are  then  deter¬ 
mined  from  (2)  by  appropriate  numerical  linear  algebra.  References 
6-8  contain  proofs  of  the  assertions  of  the  preceding  paragraph  or 
discussions  of  further  computational  considerations. 

In  contrast  to  the  situation  for  IVP's,  there  is  no  generally 
accepted  type  of  error  criterion  for  TPBVP's.  It  would  seem  desir¬ 
able  to  adopt  a  criterion  of  the  form  that  at  each  reporting  point 
at  least  one  of  the  absolute  or  relative  error  should  meet  a  speci¬ 
fied  tolerance.  The  difficulty  with  this  basic  idea  lies  in  deter¬ 
mining  a  reasonable  measure  of  the  magnitude  of  the  solution  in  de¬ 
fining  the  relative  error.  It  seems  infeasible,  or  at  least  likely 
often  to  be  inordinately  expensive,  to  take  this  magnitude  as  some 
norm  of  the  solution  at  the  particular  reporting  point  in  question. 
As  mentioned  above,  our  approach  is  to  introduce  a  class  of 
relative-error  monitors  for  the  TPBVP  (1) ,  with  each  such  monitor 
serving  to  provide  a  measure  of  the  magnitude  of  the  solution. 

More  specifically,  a  relative-error  monitor  for  (1)  is  simply 
a  real-valued  function  M=M(z)  such  that 

|  |  v(z)  |  |  <_  M(z)  .  (4) 

Given  such  a  monitor,  we  take  as  our  objective  to  control  the 
errors  <5u  and  <5v,  in  u  and  v  respectively,  so  that 

!  |  <5v  (z)  !  |  <_  e  [1+M(z)  ]  (5a) 

and 

|  |  <Su  (z)  |  |  <  £[<1+M(Z))  (1+1  |R(z)  I  |)  +  |  |qr(z)  |  |]  (5b) 

at  each  reporting  point,  where  e  (*  error  tolerance)  is  specified 
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prior  to  the  computation.  (We  shall  use  6y  to  denote  the  error 
in  a  computational  approximation  to  a  quantity  y.  We  ordinarily 
have  the  «-norm  in  mind,  although  in  the  spirit  of  our  work  it  is 
immaterial  which  of  the  common  vector  norms  is  used,  provided  the 
dependent  variables  are  properly  scaled. 

In  consonance  with  our  basic  objective,  we  wish  (5)  to  be 
implied  by  some  type  of  control  on  the  local  truncation  error 
associated  with  the  numerical  integration  of  the  IVP  (3) .  Toward 
that  end,  given  a  relative-error  monitor  for  the  TPBVP  Cl)  we  de¬ 
fine  associated  sensitivity  functions  to  be  a  quadruple  of  real¬ 
valued  functions  wR'wT»wr  and  W£  such  that  the  inequalities 

|  |  <5R (z)  |  |  <_  eWR(z)  ,  |  |  6T (z )  |  |  <_  eWT(z)  ,  (6a, b) 

I  I5qr(z)  I  I  i  eWr  (z)  ,  |  |  Sq^  (z )  j  |  <_  eW0  (z )  (6c, d) 

imply  (5).  Given  such  sensitivity  functions,  the  procedure  is  to 
integrate  the  IVP  (3)  so  that  (6)  holds  with  the  errors  interpreted 
as  the  local  truncation  errors?  if  the  IVP  is  computationally  stable 
(and  this  should  be  checked  by  multiple  numerical  integrations) , 
then  (6)  holds  with  the  errors  taken  as  the  global  errors,  and 
therefore  (5)  is  satisfied. 

All  of  the  following  are,  in  the  sense  of  the  above  definition, 
relative-error  monitors  for  (1) . 

Mq (z)  =  | | T_ 1 ( z ) | | *  1 | T (x)  |  | * | | b | |  +  |  | q^ (x) | |  +  I  I q^ (z )  |  I  ,  (7a) 

M^z)  =  i  |  T-1  (z)  T  (x)  |  ]  •  |  !  b  |  |  +  |  |T_1(z)  |  |  •  |  |q^(x)-q^(z)  |  |  ,  (7b) 

M2(z)  *  | |T_1(z)  |  | • |  |T(x)b  +  q^ (x)  — qa (*)  I  I , 


(7c) 


5. 


and 

M3(z)  *  | | v ( z ) |  j .  (7d) 

These  are  listed  in  order  of  increasing  stringency  of  the  as¬ 
sociated  requirement  on  the  error  in  v.  Sensitivity  functions 
have  been  developed  for  all  of  these  monitors.  Computational  ex¬ 
periments  currently  are  being  performed  on  a  suite  of  test  problems. 
The  results  will  be  reported  in  [9] . 

Singular  Problems 

The  method  of  invariant  imbedding  has  been  widely  applied  to 
nonsingular  TPBVP's,  such  as  (1)  of  the  preceding  subsection,  but 
little  effort  seems  to  have  been  devoted  to  its  application  to 
singular  differential  systems.  Scott  [l,Sec.V.8]  considered  suf¬ 
ficient  conditions  for  existence  of  two  different  versions  of  in¬ 
variant  imbedding  (the  "r"  and  "s"  equations)  for  linear  second- 
order  equations  having  a  regular  singular  point.  The  work  of 
Scott  constitutes  a  significant  extension  of  earlier  results  of 
Banks  and  Kurowski  [10].  More  recently,  Elder  [11]  showed  how  to 
apply  the  subscheme  of  invariant  imbedding  known  as  integration- to- 
blowup  [12,13]  to  compute  the  smallest  eigenlength  of  a  linear 
first-order  system  having  a  singularity  of  the  first  kind.  The 
purpose  of  this  work  was  to  show  how  the  approach  of  Elder  can  be 
extended  to  apply  to  the  solution  of  certain  singular  two-point 
boundary-value  problems  by  means  of  Scott's  version  [1-4]  of  in¬ 
variant  imbedding.  The  extension  to  homogeneous  differential 
systems  of  the  first  kind  has  been  completed,  as  reported  in  [14, 


15] .  Th  further  extension  to  inhomogeneous  problems  is  in  progress 
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and  will  be  described  in  [16,17].  The  application  to  homogeneous 
problems  will  now  be  described  in  some  detail. 

Consider  a  linear  homogeneous  first-order  differential 
system  with  a  singularity  of  the  first  kind  at  z=0, 


u'  (z) 
V  (z) 


+  B(z)jv(z)  , 

D0  \ 

—  +  D(z)  1  v(z)  . 


(8a) 

(8b) 


Here  u,v  and  A-D  are  as  in  the  previous  subsection,  while  Aq,Bq, 

Cq  and  Dq  are  constant  matrices  of  the  appropriate  dimensions. 

We  wish  to  obtain  solutions  of  (8)  subject  to  boundary  con¬ 
ditions  of  the  form 

u(0+)  and  v(0+)  exist  (finite)  (9a) 


and 


v(x)  =  S,  x  >  0. 


(9b) 


Let 


Y(z)  =1 


^(*)' 
kY2(z) 


(10) 


be  an  (m  +  n)  x  q  matrix  the  columns  of  which  are  a  maximal  linear¬ 
ly  independent  system  of  solutions  of  (1)  subject  to  (2a) ,  where 
Y^  and  Y2  have  respectively  m  and  n  rows.  We  define  the  (m  +  n) 


x  (m  +  n)  matrix  EQ  by 


(11) 


7. 


The  salient  results  of  [11]  can  then  be  summarized  as  follows: 

Theorem  A  (Elder  [11] J;  The  number  (q)  of  columns  in  the 
matrix  Y(z)  is  equal  to  the  dimension  of  the  space  H+  spanned  by 
the  eigenvectors  and  generalized  eigenvectors  of  E ^  associated 
with  eigenvalues  X  satisfying  Re(X)  >  0  along  with  the  ( proper) 
eigenvectors  of  E  associated  with  eigenvalue  X  =  0.  The  two- 
point  boundary -value  problem  (8),  (9)  has  a  unique  solution  if, 
and  only  if,  q  =  n  and  Y ^(x)  is  an  invertible  matrix.  If  q  =  n, 
then  Yg(x)  is  invertible  (and  hence  (8),  (9)  has  a  unique  solution) 
for  all  sufficiently  small  positive  x  if  H *  is  invertible ,  where 


is  any  (m  +  n)  x  n  matrix  the  columns  of  which  constitute  a  basis 
for  H+ .  If  there  exists  X  >  0  such  that  (1),  (2)  has  a  unique 
solution  for  every  x  satisfying  0  <  x  <  X,  then  the  matrix 


R(z)  =  Y1(z) Y 2(z)' 


satisfies  the  Riccati  differential  equation 


R  (z)  =(—■  +  B^( z)\  +  (■—  +  A^zA  R(z) 


■R(z)(AA  +  D^z^j  -  R( z)  ^ ' -j-  +  C^zl'j 


R(z) 


for  zz(0,X).  Furthermore,  if  H *  is  invertible ,  then 


8. 


R(0+)  «  H+(R+)  . 

n  n 


Henceforth  we  assume  existence  of  some  X  as  described  in  the 

statement  of  this  theorem.  The  matrices  H+  and  H+  of  (8)  are,  of 

m  n 

course,  of  respective  dimensions  m  *  n  and  n  *  n.  We  wish  to  em¬ 
phasize  that  invertibility  of  H+  implies  our  assumption  (regarding 


If  Y^  and  Y2  are  as  above,  and  u,  v  are  any  solution  of  the 

differential  system  (8)  satisfying  the  regularity  condition  (9a) , 
then  is  is  readily  shown  that  there  exists  a  constant  vector  c  such 
that 


u(z)  *  Y1(z)c,  v(z)  *  Y2U)c. 


(16a, b) 


If  additionally  v  satisfies  the  boundary-condition  (9b) ,  then  (16b) 
along  with  invertibility  of  Y2 (x)  gives 


Thus  if  we  define 


c  *  Y2(x)"1B. 


T  (z)  =  Y2(z) 


then  (16)  can  be  combined  to  give 


v(z)  *  T (z) ”XT (x)  3 ,  u(z)  =  R(z)v(z). 


(19a. b) 


9. 


Furthermore  T  is  readily  shown  to  satisfy  the  differential  system 


(z)  =  -T(z)  +  D1(z)^  +  +  C^z)^  R (z) 


(20) 


Equations  (19)  are  precisely  analogous  to  the  basic  formulas 
used  in  determining  the  solution  of  a  nonsingular  (homogeneous) 
two-point  boundary-value  problem  by  means  of  Scott's  version  of 
invariant  imbedding  [1-4] .  Furthermore  (14)  and  (15)  comprise 
an  initial-value  problem  for  R.  Thus  in  order  to  reduce  solution 
of  (8) ,  (9)  to  solution  of  initial-value  problems  and  numerical 
linear  algebra  on  matrices  and  vectors  of  order  max(m,n)  -  which 
is  the  objective  of  any  version  of  invariant  imbedding  -  it  remains 
only  to  determine  the  value  of  T(0) ,  so  that  T  can  be  obtained  by 
numerical  solution  of  the  resulting  initial-value  problem  (along 
with  that  for  R) . 

At  this  juncture  the  simple  analogy  with  the  nonsingular  sit¬ 
uation  breaks  down,  because  it  readily  can  happen  that  Y2(0)  is 
noninvertible,  so  that  T(0+)  fails  to  exist.  In  fact,  it  also 
follows  from  results  in  Ref.  11  that  Y2  (0)  is  invertible  if,  and 
only  if,  the  vector  space  H+  of  Theorem  a  consists  entirely  of 
eigenvectors  of  Eg  associated  with  the  eigenvalue  X  =  0.  In  this 
(likely  rare)  event  one  can  simply  take  T(0)  as  any  invertible 
n  x  n  matrix  (the  identity  matrix  is  certainly  the  simplest  choice) , 
because  any  such  matrix  is  Y2(0)  for  some  choice  of  Y(z).  The 
values  of  R  and  T  can  then  be  obtained  by  numerical  integration  of 
the  initial-value  problem  consisting  of  (14)  ,  (15)  and  (20)  ,  sub¬ 
ject  to  the  initial  value  selected  for  T. 
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In  the  case  that  Y2(0)  is  singular,  the  only  general  procedure 
we  have  found  is  as  follows.  It  is  readily  shown  that  satisfies 
the  differential  equation 


Y' (z) 


C1  (z) 


R(z) 


+ 


Dl(z) 


) 


Y2(z)  . 


(21) 


For  0  <_  z  <_  Zq  one  can  then  determine  R(z)  and  Y2(z)  by  numerical 
integration  of  the  initial-value  problem  consisting  of  (.14)  ,  (15)  , 
(21)  and  the  initial  condition 


Y2(°>  =  P2  ' 


(22) 


where 


is  any  (m  +  n)  *  n  matrix  the  nonzero  columns  of  which  constitute 
a  basis  for  the  null  space  of  Efl.  With  Y2  (zQ)  thus  determined, 
T(zq)  =  Y2(zq)  1  can  be  obtained  numerically.  The  required  values 
of  R(z)  and  T(z)  for  z  >  Zq  can  then  be  obtained  by  numerical  in¬ 
tegration  of  (14)  and  (21)  subject  to  the  known  values  of  R(zq) 
and  T(zq).  Alternately  one  could  simply  compute  Y2(z)  and  R(z) 
for  all  desired  values  of  z,  and  then  compute  u  and  v  by 


v  ( z )  =  Y2(z)Y2(x)_1B 


(23) 


and  (19b) . 
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There  remains  one  further  computational  detail.  Typical 
initial-value  codes  require  the  user  to  supply  a  subroutine  which 
computes  the  derivative.  However,  the  right-hand  side  of  (14)  is 
indeterminate  at  z  *  0  for  R(0)  given  by  (15)  and  (12) ,  and  like¬ 
wise  the  right-hand  of  (20)  (respectively  (21))  in  the  case  that 
T(0)  (respectively,  ^  (0) )  exists.  For  the  problems  we  have  en¬ 
countered  it  has  been  possible  to  determine  R'  (0)  and  Y^(0) 

(or  T ' ( 0 ) ,  when  T(0)  exists)  by  applying  L'Hospital's  rule  to  (14) 
and  (21)  (or  (20)),  although  we  have  no  general  proof  of  the  effec¬ 
tiveness  of  this  approach. 

Several  numerical  examples  illustrating  the  approach  outlined 
above  are  presented  in  references  15  and  16.  In  addition  the  im¬ 
portant  special  case  of  a  single  second-order  equation  is  further 
studied  in  considerable  detail,  and  several  results  specific  to 
such  problems  are  obtained. 
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